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Nomenclature 

normal  boundary  layer  coordinate 
streamwise  boundary  layer  coordinate 
dimensionless  stream  function  variable 
specified  wall  shear 

number  of  steps  in  the  normal  direction 
number  of  steps  in  the  streamwise  direction 

right  hand  side  groups  for  spline  correction 
equations  containing  known  quantities  from 
previous  iteration  or  streamwise  stations 


4x4  block  tridiagonal  coefficient  matrices 

column  vector  of  known  quantities 
solution  correction  vector 
convergence  criteria  value 
iteration  number 
streamwise  grid  location 
normal  grid  location 


All  other  quantities  are  defined  in  the  text. 
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I.  Introduction 

The  object  of  the  present  study  is  the  accurate  numerical  determination  of 
the  pressure  gradient  distribution  for  specified  wall  shear  in  a  laminar 
boundary  layer.  This  is  known  as  the  inverse  problem.  This  procedure  is 
sometimes  preferable  to  the  standard  problem,  in  which  the  pressure 
distribution  is  prescribed  and  the  body  shear  stress  is  determined  by  the 
solution  of  the  boundary  layer  equations. 

One  method  of  determining  the  pressure  gradient  distribution  is  the 
mechul  function  scheme,  developed  by  Cebeci  and  Keller  [1].  This  scheme 
was  found  to  be  accurate  and  efficient  for  self-similar  type  flow  problems. 

For  non-similar  flow  problems,  however,  this  method  exhibits  a  weak 
numerical  instability.  The  instability  appears  as  a  disturbance  that 
develops  in  the  far  field  of  the  computational  domain  and  propagates  towards 


the  wall.  As  a  result,  the  numerical  solution  is  ultimately  destroyed  as 
It  is  marched  downstream. 

The  Inverse  problem  requires  the  determination  of  a  coefficient  function, 
the  pressure  gradient  parameter,  which  satisfies  an  over-determined  set  of 
boundary  conditions  that  result  when  the  wall  shear  is  specified.  The  mechul 
function  scheme  adds  an  extra  differential  equation  for  the  pressure  gradient 
to  the  system  of  boundary  layer  equations.  As  a  result,  the  system  is  no 
longer  over-determined,  and  the  governing  equations  may  be  solved  as  an 
extention  of  the  standard  problem. 
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This  study  reformulates  the  mechul  function  scheme  as  developed  in  Ref. 

[1]  in  such  a  way  that  the  numerical  instability  is  eliminated,  thus  yielding 
an  accurate  solution  method.  Whereas  in  Ref.  [1],  the  Keller  box  scheme  was 
used  to  discretize  the  differential  equations,  the  new  procedure  applies 
fourth  order  splines  developed  by  Rubin  and  Khosla  [3]  to  approximate  the 
normal  derivatives  and  three  point  backward  differences  for  the  streamwise 
derivatives.  This  yields  a  fully  implicit  method  which  lends  stability  to  the 
solution. 

In  the  present  work,  the  discretized,  linearized  system  of  equations 
put  into  a  standard  block  tridiagonal  form.  Lower-upper  decomposition  .  ’d 
to  solve  the  resulting  block  matrix  system  with  partial  pivoting  within  toe 
blocks  to  prevent  the  buildup  of  round-off  errors.  Reformulation  of  the 
governing  equations  became  necessary  with  the  discovery  that  the  spline 
techniques,  applied  to  the  original  formulation  in  Ref.  [1],  yielded  a 
singular  block  matrix  on  the  diagonal  during  forward  substitution.  By  a 
slight  recasting  of  the  differential  equation  governing  the  pressure  gradient, 
this  singularity  is  removed  in  the  matrix. 

This  study  considers  the  solution  of  both  self-similar  and  non-similar 
laminar  boundary  layer  problems.  Falkner-Skan  problems  with  positive  and 
negative  wall  shear  and  two  non-similar  flow  situations  are  computed. 

II.  Analysis 

Governing  Differential  Equations 

The  governing  equations  in  this  problem  are  the  two-dimensional  boundary 
layer  equations  for  incompressible  laminar  flow.  In  dimensionless  stream 
function  form,  the  boundary  layer  equation  is  written 
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Equations  (L),  (3)  and  (4)  lead  to  an  over-determined  system.  To  complete 
the  formulation  of  the  raechul  function  scheme,  the  pressure  gradient  parameter 


is  written  as 


0(0  =  8(s,n)  . 


Then  in  Ref.  [1]  the  following  derivative  condition  is  introduced: 


M-o 


K  >  0  . 


The  pressure  gradient  parameter  is  thus  determined  by  solving  Eqs.  (1)  and 
(5)  together  with  boundary  conditions  (3)  and  (4). 

For  self-similar  flows, 

f5(C,n)  =  0  .  (6) 

Therefore,  the  boundary  conditions  and  the  pressure  gradient  are  independent 
of  the  streamwise  coordinate.  In  this  case,  the  differential  equations 
reduce  to 


f  +  ff  +  g(l  -  f  1  =  0 
rmn  nn  v  n' 


8=0  (7b) 

n 

yielding  a  system  of  ordinary  differential  equations. 

It  was  found  during  the  present  study  that  the  original  gradient  parameter 
condition,  Eq.  (5),  used  with  the  fourth  order  spline  relations,  leads  to  a 
singular  matrix  when  solving  the  resulting  block  tridiagonal  system  by 
L-U  decomposition.  Using  a  constant  normal  stepsize  in  the  spline  S^(4,0) 
given  in  Ref.  [3]  yields  a  singular  block  matrix  on  the  diagonal,  thus 
preventing  the  solution  of  the  matrix.  This  problem  became  apparent  in 
the  forward  substitution  step. 
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To  eliminate  this  singularity,  the  pressure  gradient  parameter  condition 
was  changed  to 


328 


=  0 


£  >  0 


(8) 


3rf 


Then,  the  condition 


5  >  0 


(5) 


is  enforced  at  the  wall  with  two  point  spline  boundary  conditions.  Using  the 
alternate  fourth  order  spline  relation  S^(4,0)  from  Ref.  [3]  produces  a  non¬ 
singular  matrix  for  constant  or  varying  stepsize  while  still  forcing  the 
pressure  gradient  to  be  independent  of  n. 


Differential  Equations  in  First  Order  Form 

The  governing  differential  equations,  (1),  are  written  in  first  order  form 
with  the  exception  of  Eq.  (8): 


f  =  u 

n 


u  =  T 
n 


8  =  0 
nn 


(9a) 

(9b) 

(9c) 


T  =  8(l  -  u  )  -  ft  +  2?(uu-  -  Tf _) 


(9d) 


Then  the  boundary  conditions  are 


f(c,0)  =  0 

(10a) 

u(5,0)  =  0 

(10b) 

i 

g 
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t(C,0)  =  S(C)  (10c) 


u(c,nj  =  l  .  (lOd) 


Streamwise  Discretization 

The  governing  differential  equations  are  discretized  using  fourth  order 
accurate  splines  in  the  normal  direction  and  backward  differences  in  the 
streamwise  direction.  For  the  £  derivatives,  the  following  general  backward 
difference  formula  for  constant  stepsize  is  used. 

For  i  >  2,  the  three  point  second-order  accurate  version  is  used  with 


a  =  |  ,  b  =  -  1  ,  and  c  -  ~ 


(11a) 


whereas,  for  the  second  downstream  £  station  (i  =  2),  the  first-order  accurate 


(two  point)  formula  is  used  with 


a  =  1  ,  b  =  -  1  ,  and  c  =  0 


(lib) 


Examining  the  governing  equations,  only  Eq.  (9d)  contains  streamwise 
derivaties.  Applying  Eq.  (11)  to  the  and  f^  derivatives  yields 

=  Bi,jK,j " _  fi,.iTi,j 


+  °i,J(auLj  +  bui-l,3  +  CVl,j] 


‘  2“i,j(afi,jTi,j  +  bfi-l,jTi,j  +  Cf  i-2 ,  j  Ti  ,  j  ^ 


Tarawa  mmt  mb  ^ 


rV  ■  '  .*  V"  •  •  •.'■  •  •  ’  -  V  ■ 
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where 


t,j  4C4 


Normal  Discretization 

Since  first  and  second  derivatives  with  respect  to  n  occur  in  the 
equations,  the  splines  S^(4,0)  and  S^(4,0)  are  used.  To  apply  the  spline 
approximations  in  the  normal  direction,  the  following  spline  derivatives 
are  defined: 


Bnn  * 


These  definitions  are  then  substituted  into  the  governing  equations, 
Substituting  definitions  (13)  into  Eqs.  (9)  yields 


*1,3  "  Ui,3 


»u  =  r 

i.j  i,J 


(13a) 


(13b) 


(13c) 


(13d) 


(I3e) 


(14a) 


(14b) 
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( 14c) 


+  +  +  CU?-2,.,)  • 


(14d) 


From  Ref.  [3],  the  tridiagonal  relation  for  Sl(4,0)  is,  at  a  streamwise 


station  i. 


<5+1  Mi  ♦  .)2««  ♦  A*  ■ 


2  1  r  1  +  2  a 

1  +  a  hj  l  0 


t  +  slz± 

’j+l  o 


+  o)3gi 


-  o  (2  +  ojgj.J 


hj '  4nj  • 


, . „ . 

j  hj  ’ 


where 
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and  is  the  first  derivative  spline  approximation  of  (3g/3n)j*  The 
tridiagonal  relation  for  S^(4,0)  is  given  as 

2  3  2  2 

<J  ~  a  ~  1  .g  ,  o  +  4q  4q  -t-  1  ,g  ,  l+o-o  Tg 

12a  j+1  12a  j  12  Lj-1 


=  L.  [fill  -  l+_2  g  +  g  1 

.  2  L  a  a  8 j  gj-lJ 


(16) 


j 


.2  ,2 


where  Lj  is  the  second  derivative  spline  approximation  of  (3  gj  3n  )j • 

To  eliminate  the  spline  derivatives,  Eqs.  (14a),  (14b)  and  (14d)  are 
substituted  into  the  S^(4,0)  relation  given  by  Eq.  (15).  The  second 
derivative  condition,  Eq.  (14c),  can  be  substituted  into  the  S^(4,0) 
relation.  This  yields  the  following  four  tridiagonal  relations 
(i  subscript  understood).  All  coefficients  are  given  in  Appendix  A. 


o2u.  ,  +  alu,  +  u, ..  =  o3f,  ,  +  o4f,  +  o5f... 

3-1  i  J+1  j"1  1  J+1 


o2x,  .  +  air,  +  x...  =  a3u.  .  +  a4u.  +  a5u. 
j-1  J  J+1  J“1  J  J+1 


°86j_1  +  o7  8j  +  a6  Bj =  0 


(17a) 


(17b) 


(17c) 


W.  .It  ■  >-3-^. 
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o2 [ (cl  +  8)u2  -  8  +  c2f t  +  El. t' 


j-1 


+  ol[(cl  +  B)u  -  8  +  c2fx  +  E2^.  t] 


.1 


0  _ 

[(cl  +  B)u  -  8  +  c2f t  +  E3 . x] 


J  j+1 


=  -  o2C4.  ,  -  olc4  -  c4 

j-1  j  J+l 


Linearization 

The  block  tridiagonal  relations  given  by  Eq.  (17)  form  a  nonlinear 
system  relating  (f ,u, t.B)^ j •  These  equations  are  first  linearized  and 
then  solved  using  Newton's  method.  The  Newton  iterates  are  given  by 


f(n+l) 


+  6f 


(n) 


( 


(n+1) 

Ui,j 


+  6u 


(n) 


( 


i 


T 


(n+1) 


(n 

i. 


+  6t 


(n) 


( 


i 

i 

i 

i 

i 

i 


B 


(n+1) 


,(n) 

l,j 


+  58 


(n) 


( 


where  the  superscript  indicates  the  iteration  number.  Equations  (18)  are 
substituted  into  the  tridiagonal  relations  at  the  (n+l)st  iteration  and 
the  quadratic  and  higher  order  terras  are  neglected.  This  allows  the 


A'9  ta\ t'M'S* 
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equations  to  be  solved  for  the  unknown  corrections,  with  all  other  terms 
being  known  at  iteration  (n).  Moving  known  quantities  to  the  right  hand 
side,  the  linearized  correction  equations  in  block  tridiagonal  form  at 
streamwise  station  i  are 


-  o36f .  .  +  o26u 


.  .  -  a46f.  +  ol6u.  -  o56f.,.  +  6u...  =  L.  (19a) 

J“l  J  J  j+1  j+1  J 


o36uj_^  +  o26tj_^  -  o46uj  +  oIStj  -  o56uj  +  1  +  6tj  +  1  =  Pj  (19b) 


3j_1  -  a7  60j  -  a66Bj+1  =  Qj 


(19c) 


o2[A6f  +  B6u  +  D^St  +  G5e]j_1  +  ol  [A6f  +  B6u  +  E^  St  +  G6b]j 

A  A  A  A 

+  [A6f  +  BSu  +  Fj6t  +  G«b]j+1  “  Tj 


(19d) 


where  A j ,  B  j ,  D j ,  Ej  and  Fj  are  coefficients  obtained  from  the  linearization 
of  the  equations  (refer  to  Appendix  A). 

These  block  tridiagonal  equations  can  then  be  written  in  the  following 
matrix  form  (i  subscript  understood) 

Vj-1  +  Vj  +  CjZj+l  "  Rj  »  2  *  j  <  N  (20) 

where 


6f 

6u 

St 

SB 

- 


\  ■  ,  y„r  »r,  v\  .*  ,  *”*.•.  •  .  ".  "  .  *.  .*.  .  *  -* 
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and  A j ,  B j ,  Cj  are  4x4  block,  matrices  and  Rj  is  the  column  vector  of 

quantities . 


known 


Two  Point  Spline  Relations  at  the  Wall 


The  boundary  conditions  at  the  wall  in  correction  form,  at  streamwise 
station  i,  are 


6f  L  =  0 


6u^  =  0 


6t^  -  0 


One  additional  condition  must  be  provided  to  close  the  system  at  the  wall, 
two  point  spline  boundary  condition,  given  in  Ref.  [3],  is  used. 


b2  -  »1  -  r  (»2*  -  »?) +  t§  (hs  -  hB)  - 0  • 


This  relation  allows  the  original  pressure  gradient  condition  to  be  enforced 
at  the  wall.  From  the  governing  equations, 

L?  -  0 


and,  enforcing  the  original  gradient  condition 


‘j-0  ■ 


yields  the  simple  relationship 


82  =  81 


or,  in  correction  form 
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682  -  <S81 


(24b) 


The  linearized  block  system  for  the  solution  corrections  at  the  wall  can  then 
be  written  in  the  following  matrix  form: 


A1Z1  +  C1Z2  =  R1 


(25) 


with  Zj  defined  by  Eq.  (21). 

Two  Point  Spline  Relation  at  the  Outer  Edge 
The  far  field  boundary  condition  is,  at  any  i, 


6uN+l  “  0 


(26) 


Three  additional  conditions  are  required  to  close  the  system.  The  two  point 
spline  conditions  are  again  applied.  For  one  relation,  the  following  is 
used: 


N+l  N 


,  hN+l  (  .  >,  .  hN+l  ( 

f”  “  9  (UN+1  +  UN'  +  19 


12  '‘‘N+l  TN^  0 


(27a) 


Differentiating  Eq.  (27a)  with  respect  to  n  yields 


N+l 


L  I.*1 

N+l  ,  ^  |  N+l  r„  T  _ 

"  UN  “  2  ^TN+1  +  V  +  12  ^N+i  0 


(27b) 


Differentiating  once  more  gives 


N+l 


hN+l  (  ^  . t ^  hN+l  T 

“  T”  9  v  ^N+l  '  19  (L*J4 


N 


12  v  N+l 


-  O  "  0 
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Assuming  the  difference  between  the  second  order  terms  to  be  negligible 
compared  to  the  lower  order  terras,  the  above  equation  simplifies  to 


T  —  T  — 

N+l  N 


K+l  +  *nT)  ■  0  • 


(27c) 


Equations  (27)  provide  the  three  additional  conditions  required  for  the  far 
field.  These  equations  are  linearized  and  solved  for  the  correction  terras  as 
before  and  then  written  in  matrix  form. 


VlZN  +  AN+1ZN+1  Vt-l 


Solution  of  Block  Tridiagonal  Matrix  Equation 

The  block  tridiagonal  system  formed  by  Eqs.  (20),  (25)  and  (28)  is  solved 
using  standard  lower-upper  (LU)  decomposition  together  with  partial  pivoting 
within  the  4x4  blocks  to  prevent  the  buildup  of  roundoff  errors.  A  block 
tridiagonal  solver  using  subroutines  developed  by  Blottner  [4]  which  perform 
the  partial  pivoting  is  used  to  solve  the  matrix  equations  for  the  corrections 
at  each  iteration. 


Starting  Solution 


The  Falkner-Skan  self-similar  equations  are  obtained  by  setting  £  equal  to 


zero  in  Eqs.  (1)  and  (8). 


f  +  ff  +  Sfl  -  f  1  =  0 

nnn  nn  v  n' 


0  =  0  . 
nn 


(29a) 


(29b) 


The  boundary  conditions  given  by  Eqs.  (3)  and  (4),  with  no  dependence  on 
are  used.  The  solution  of  this  ordinary  differential  equation  with  splines 


* .  •  <1  •  •  .<  . 
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~  1 

Isi 

l 


approximating  all  n  derivatives  is  used  as  the  starting  solution  for  the 
non-similar  case.  The  scheme  is  derived  so  that  the  Falkner-Skan  solutions 
for  positive  and  negative  wall  shear  can  also  be  computed. 

4s  an  initial  guess  for  the  starting  solution,  a  fourth  order  Pohlhausen- 
type  polynomial  is  used  to  approximate  ( f ,u, t, 8)t f j  at  5=0.  The  Pohlhausen 
polynomial  for  the  velocity  is  of  the  form 

u  =  +  c^  +  d£^  +  e^  (30) 


where 


?  =  -  ,  u  =  u(  n)  , 


and  the  constants  b,  c,  d  and  e  can  be  found  by  applying  the  boundary 
conditions  and  the  ordinary  differential  equation.  The  stream  function  and 
the  shear  can  be  found  by  integrating  and  differentiating  Eq.  (30) 
respectively.  By  substituting  into  the  differential  equation,  Eq.  (29a), 
an  approximation  for  8  is  obtained  in  the  form 


6Sn  -  12 


Once  the  starting  solution  is  determined,  the  second  streamwise  station 
must  be  treated  in  a  special  way  for  the  non-similar  flow  case.  With  only 
one  previous  streamwise  step  known,  two  point  backward  differences  are  used 
to  approximate  the  ?  derivatives.  Past  this  station,  three  point  backward 
differences  are  used.  The  solution  profile  at  the  last  calculated  station 
is  used  as  the  first  approximation  at  the  new  station. 


wcwe 


V.V 


to 


59 
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II .  Results  and  Discussion 

Similar  Flow  Problem 

Computations  for  similar  flows  were  performed  for  positive  and  negative 
wall  shears.  For  all  positive  wall  shears,  the  solutions  were  obtained 
independently  for  a  specific  shear.  For  the  negative  wall  shears,  the 
sensitivity  of  the  method  to  the  initial  guess  required  calculating  solutions 
consecutively  for  small  steps  in  shear  and  using  the  previous  solution  profile 
as  an  initial  guess  for  the  next  profile.  The  solution  process  was  begun  for 
a  zero  wall  shear,  S  =  0,  and  the  shear  was  decremented  by  0.01  or  0.05. 

Solution  comparisons  are  made  between  the  present  reformulated  scheme, 
the  original  mechul  function  formulation  [1],  and  the  nonlinear  eigenvalue 
scheme  [2].  This  method  was  developed  by  Keller  and  Cebeci  before  the  mechul 
function  scheme.  The  eigenvalue  method  solves  the  inverse  problem  by  treating 
the  unknown  pressure  gradient  as  an  eigenvalue.  Then,  two  iteration 
procedures,  an  "inner”  and  an  "outer"  iteration,  are  performed.  The  inner 
iteration  solves  the  governing  equations  for  a  standard  problem  assuming  8 
is  known.  This  inner  iteration  is  then  used  with  Newton's  method  to  determine 
the  pressure  gradient  parameter  using  the  variational  equations  in  the  outer 
iteration  procedure.  The  variational  equations  are  the  standard  boundary 
layer  equations  differentiated  with  respect  to  B. 

For  positive  wall  shears  computed  with  the  reformulated  scheme,  a  normal 
stepsize  of  An  =  0.15  is  used  with  Pco  =  6;  for  reverse  flow  solutions  and 
separation  (S  =  0),  An  =  0.15  and  n*  =  9  are  used.  The  criteria  used  for 
convergence  is 

I  .(n+ll  (nil  ..-8 
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Cebeci  and  KelLer  appLled  a  similar  convergence  test  to  the  original  scheme  as 
well  as  the  eigenvalue  scheme  with  e  <  lO-^  [1,2]. 

The  results  of  the  self-similar  flow  calculations  with  positive  wall 
shear  are  given  in  Table  I.  These  results  are  compared  with  those  of  Smith 
[5].  Comparison  shows  the  values  from  the  reformulated  mechul  function  scheme 
closely  approach  those  of  Smith.  All  cases  converged  quadratically .  It 
should  be  noted  that  the  greater  number  of  iterations  required  for  convergence 
with  the  reformulated  scheme  are  the  result  of  the  present  more  severe 
convergence  criteria. 

The  results  of  the  reverse  flow  computations  are  presented  in  Table  II. 
These  results  are  compared  with  those  of  Stewartson  [6],  Again,  the  agreement 
is  good.  Here,  the  number  of  iterations  decreases  appreciably  with  the  use  of 
consecutive  calculations.  Figure  2  shows  examples  of  both  positive  and 
negative  wall  shear  velocity  profiles. 

Non-Similar  Flow  Problem 

Non-sirailar  flow  computations  were  performed  using  two  linear  wall  shear 
distributions.  The  first  case,  given  by 

S(C)  =  0.4696  (1  -  C) 

has  a  zero  pressure  gradient  at  5=0,  indicating  a  flat  plate  flow.  The 
second  case,  given  by 

S(C)  =  1.23259  (  1  -  C) 

has  a  pressure  gradient  parameter  of  unity  at  5=0,  indicating  a  stagnation 
point.  Both  cases  approach  zero  shear  at  5  =  1,  yielding  separation.  Both 
cases  were  computed  with  values  of  n«,  =  6  and  n®  =  9.  For  all  non-similar 
computations,  a  strearawise  stepsize  of  AC  =  0.05  was  used,  with  An  =  0.25. 
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The  convergence  criteria  applied  is  identical  to  that  used  for  the  self- 
similar  cases. 

A  comparison  of  results  for  the  first  case  is  given  in  Tables  III  and  IV. 
Table  III  compares  the  results  with  those  of  Refs.  [1]  and  [2].  The 
agreement  with  the  Richardson  extrapolation  results  obtained  from  the 
eigenvalue  scheme  values  is  quite  good.  The  numerical  instability  experienced 
by  the  original  formulation  does  not  appear  in  the  present  method.  It  was 
found  that  with  n<»  =  6 ,  the  solution  would  not  converge  at  £  stations  near 
separation,  yet  with  n<»  =  9,  the  solution  would  march  through  the  separation 
point  at  5  =  1.  Once  past  separation,  however,  the  scheme  quickly  becomes 
unstable  and  the  solution  does  not  converge.  Table  IV  compares  results  for 
Case  1  with  n®  =  6  and  n®  =  9. 

A  comparison  of  results  for  the  second  case  is  given  in  Tables  V  and  VI. 

In  Table  V,  the  results  are  compared  with  results  from  Ref.  [2],  since  no 
results  are  available  for  this  case  with  the  original  mechul  function 
formulation.  The  agreement  with  the  Richardson  extrapolation  is  again  very 
good.  As  with  case  one,  the  results  could  be  obtained  through  separation  only 
with  n®  =  9 .  Table  VI  compares  the  results  for  the  two  values  of  nM. 

Figures  3  and  4  are  plots  the  pressure  gradient  parameter  as  a  function  of 
€  for  both  non-similar  cases.  The  results  are  smooth,  with  a  minimum 
occurring  just  before  separation  for  computations  performed  with  =  9. 

A  Note  on  Normal  Stepsize 

For  all  numerical  cases  presented,  the  normal  stepsize,  An,  is  kept 
constant  for  reasons  of  comparison  with  Refs.  [1]  and  [2].  Non-similar 
cases  were  computed,  however,  using  a  geometric  progression  for  n,  with 
hj+l/hj  3  1.1.  Using  this  geometric  progression,  both  non-similar  cases 


-24-  4  May  1984 

KCK:GHH: Ihz 

would  proceed  through  separation  with  n®  =  6 .  The  results  obtained  are 
identical  with  those  produced  with  a  constant  stepsize  and  n®  =  9.  Using 
the  geometric  progression  also  allowed  the  number  of  normal  steps  to  be 
halved  without  decreasing  solution  accuracy. 

IV.  Conclusion 

An  accurate,  stable,  and  efficient  method  is  developed  to  determine  the 
pressure  gradient  distribution  on  a  body  surface  in  a  laminar  boundary  layer 
flow  with  wall  shear  specified.  The  reformulated  mechul  function  scheme 
presented  here  does  not  exhibit  the  numerical  instability  experienced  with 
non-similar  type  flow  problems  in  Ref.  [1].  The  method  is  fully  implicit 
and  is  applicable  to  Falkner-Skan  as  well  as  non-similar  flow  problems. 

Reverse  flow  self-similar  problems  have  also  been  computed  but  are  found 
to  be  very  sensitive  to  the  initial  solution  guess. 

The  reformulated  scheme  uses  fourth  order  splines  to  approximate 
derivatives  in  the  normal  direction  and  three  point  backward  differences  for 
the  streamwise  derivatives,  both  to  aid  solution  stability.  Partial  pivoting 
is  used  within  the  4x4  blocks  resulting  from  the  discretized,  linearized 
equations  of  motion.  These  modifications  prevent  the  instability  described  in 
Ref.  [1]. 

The  results  for  both  similar  and  non-similar  cases  were  found  to  be  stable 
and  accurate  with  quadratic  convergence  consistently  observed.  It  was  found 
that  solution  calculations  would  proceed  through  the  separation  point  with 
constant  stepsize  and  n*  =  9,  or  with  a  geometric  progression  and  n®  =*  6. 
Although  all  solutions  converged  quadratically  upstream  of  the  separation 
point,  using  a  more  accurate  initial  guess  would  decrease  iteration  counts 


further. 
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TABLE  3E 

COMPUTED  PRESSURE  GRADIENT  MS  A  FUNCTION  OF  £  FOR  THE  WALL 
SHEAR  DISTRIBUTION: 

S(£)  =  1.23259  (K) 
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Figure 


Computational  Domain 


PRESSURE  GRADIENT  PARAMETER  AS  A  FUNCTION  OF  £ 


Figure  4 

PRESSURE  GRADIENT  AS  A  FUNCTION  OF  S 
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Appendix  A 


Fourth  Order  Spline  Coefficients: 
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Streamwise  Discretization  Coefficients; 
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